Decomposotion Methods

Benders Decomposition

Suppose we have an optimization problem given by

\begin{align} \tag{P} \begin{array}{ll} \min & c^\intercal x + f(y)\\ \text{s.t.} & Ax + F(y) \geq b \\ & x \geq 0, \ y \in Y, \end{array} \end{align}

where A is m \times n, f:\Re^q \mapsto \Re, F:\Re^q \mapsto \Re^m and Y \subseteq \Re^q. The difficulty in this problem comes from the variable y as the problem becomes a linear program for a fixed y. In a typical application the set Y can also be complicated set like the integrality constraints. This type of formulation also arises frequently in stochastic programming, where the decision variables y are the first stage variables, and variables x represent the activities in the future. As the uncertainty in the future reveals itself, the variable y is altered and the corresponding expected cost, denoted as f(y), is minimized.

Let \mathscr{Y} = \{y \in Y: \exists x \ge 0 \text{ with } Ax \ge b - F(y)\}, and \mathscr{U} = \{A^\intercal u \le c, u \ge 0\}. We rewrite (P) as follows:

\begin{align*} &\min_{x, y \in Y} \bigg\{ f(y) + c^\intercal x : Ax \ge b - F(y), x \ge 0 \bigg\} \\ \stackrel{[1]}{=} &\min_{y \in \color{red}{\mathscr{Y}}} \bigg\{ f(y) + \min_x \left\{ c^\intercal x: Ax \ge b - F(y), x \ge 0 \right\} \bigg\} \\ \stackrel{[2]}{=} &\min_{y \in {\mathscr{Y}}} \bigg\{ f(y) + \max_u \left\{ (b - F(y))^\intercal u: A^\intercal u \le c, u \ge 0 \right\} \bigg\} \\ \stackrel{[3]}{=} &\min_{y \in \mathscr{Y}} \bigg\{ f(y) + \max_{i \in K} \left\{ (b - F(y))^\intercal u^p_i \right\} \bigg\} \\ = & \min \quad f(y) + z \\ & \begin{aligned} \text{ s.t.} \quad z &\ge (b - F(y))^\intercal u^p_i, && i \in K \\ y & \in \mathscr{Y} \end{aligned} \\ \stackrel{[4]}{=} & \min \quad f(y) + z \\ & \begin{aligned} \text{ s.t.} \quad & z \ge (b - F(y))^\intercal u^p_i, && i \in K \\ & \{x \in \mathcal{R}^n_{+} : Ax \ge b - F(y)\} \neq \emptyset && \end{aligned} \\ \stackrel{[5]}{=} & \min \quad f(y) + z \\ & \begin{aligned} \text{ s.t.} \quad z &\ge (b - F(y))^\intercal u^p_i, && i \in K \\ 0 & \ge (b - F(y))^\intercal v^r_i, && i \in R \\ y & \in Y \end{aligned} \end{align*}

• The equation [1] is simply a rewrite of (P). Note that y \in \mathscr{Y} guarantees the feasibility of the LP on x.
• The equation [2] holds because of the LP’s duality.
• The equation [3] holds because of u^p_i, i \in K, are extreme points of the polyhedral \mathscr{U}.
• The equation [4] follows the definition of \mathscr{Y}, which indicates that y \in \mathscr{Y} is equivalent to \{x \in \mathcal{R}^n_{+} : Ax \ge b - F(y)\} \neq \emptyset.
• The equation [5] holds because by the Farkas’ Lemma, \begin{align*} \{x \in \mathcal{R}^n_{+} : Ax \ge b - F(y)\} \neq \emptyset \Leftrightarrow \ & \text{$(b-F(y))^\intercal v \le 0$ for all $v \in \Re_+^m $ that $ A^\intercal v \le 0$} \\ \Leftrightarrow \ & \text{$(b-F(y))^\intercal v \le 0$ for all extreme ray $v$ of the polyhedral cone $\{v \in \Re_+^m : A^\intercal v \le 0\}$} \\ \Leftrightarrow \ & \text{$(b-F(y))^\intercal v^r_i \le 0$ where $v^r_i$, $i \in R$, are extreme rays of $\mathscr{U}$} \end{align*}

Capacitated Fixed-Charge

We consider a capacitated fixed-charge problem with 4 suppliers and 8 customers as shown in the figure

The mathematical model is:

\begin{align*} \min \ & \sum_{j \in S} f_j y_j + \sum_{i \in C, j \in S} t_{ij} x_{ij} \\ \text{s.t.}\ & \sum_{j \in S} x_{ij} = d_i, \ \forall i \in C, \\ & \sum_{i \in C} x_{ij} \le c_j y_j, \ \forall j \in S, \\ & ~~ x_{ij} \ge 0,\ y_j \in \{0,1\} \end{align*}

We have

Master Problem (Initial)	Subproblem (Primal)	Subproblem (Dual)
\begin{align*} \min \ & \sum_{j \in S} f_j y_j + z \\ \text{s.t.}\ & \sum_{j \in S} c_j y_j \ge \sum_{i \in C} d_i \\ & ~~ \ y_j \in \{0,1\} \end{align*}	\begin{align*} \min \ & \sum_{i \in C, j \in S} t_{ij} x_{ij} \\ \text{s.t.}\ & \sum_{j \in S} x_{ij} = d_i, \ \forall i \in C \\ & \sum_{i \in C} x_{ij} \le c_j y_j, \ \forall j \in S \\ & ~~ x_{ij} \ge 0, \end{align*}	\begin{align*} \max \ & \sum_{i \in C} d_{i} u_{i} - \sum_{j \in S} c_j y_j v_j \\ \text{s.t.}\ & u_i - v_j \le t_{ij},\ \forall i \in C, j \in S \\ & v_j \ge 0 \end{align*}

and the restricted Benders master problem \begin{align*} \min \ & \sum_{j \in S} f_j y_j + z \\ \text{s.t.}\ & \sum_{j \in S} c_j y_j \ge \sum_{i \in C} d_i \\ & \sum_{i \in C} d_{i} u^i_{i} - \sum_{j \in S} c_j y_j v^i_j \le z,\ \forall i \in \bar{K} && \text{optimality cuts} \\ & \sum_{i \in C} d_{i} u^r_{i} - \sum_{j \in S} c_j y_j v^r_j \le 0,\ \forall r \in \bar{R} && \text{feasibility cuts} \\ & \ y_j \in \{0,1\} \end{align*}

Note

The constraint in the master problem, used capacity greater or equal to total demand, is not necessary, but it helps to tight the LP relaxation.

Without this constraint, we will get y_j = 0 as the optimal solution in the first iteration. This solution will cause the subproblem infeasible (or unbounded for the dual problem), and feasibility cut is required to be added into the master problem.

The direct formulation of the dual is

\begin{align*} \max \ & \sum_{i \in C} d_{i} u_{i} + \sum_{j \in S} c_j y_j v_j \\ \text{s.t.}\ & u_i + v_j \le t_{ij}, \ \forall i \in C, j \in S \\ & v_j \le 0 \end{align*}

This is the dual of the primal problem used in Gurobi. We formulate the dual problem by substitute -v_j as v_j. Hence, we need to add negative sign in the implementation after we solve the primal and obtain the dual solution.

Code

import pandas as pd
import seaborn as sns
import numpy as np
import matplotlib.pyplot as plt
from gurobipy import *

sns.set_style('ticks')
%matplotlib inline
plt.rcParams['figure.dpi'] = 70

fixed_cost = {'S1':4500, 'S2':4400, 'S3':4250, 'S4':4250}
capacity   = {'S1':600,  'S2':700,  'S3':500,  'S4':550}
demand     = {'D1':100,  'D2':150,  'D3':175,  'D4':125, 
              'D5':180,  'D6':140,  'D7':120,  'D8':160}
trans_cost = [26, 27, 28, 17, 13, 19, 21, 19, 19, 22, 12, 
              26, 11, 22, 22, 24, 17, 15, 13, 15, 25, 26, 
              27, 23, 20, 19, 23, 26, 21, 16, 23, 26]
suppliers = list(fixed_cost.keys())
customers = list(demand.keys())
trans_cost = dict(zip([(i,j) for i in customers for j in suppliers],trans_cost))
# np.array(list(trans_cost.values())).reshape(len(customers), -1).T

m = Model()
y = m.addVars(suppliers, vtype=GRB.BINARY, name='y')
x = m.addVars(customers, suppliers, vtype=GRB.CONTINUOUS, name='x')

m.setObjective(sum(fixed_cost[j]*y[j] for j in suppliers) + 
               sum(trans_cost[i,j]*x[i,j] for i in customers for j in suppliers))

for i in customers:
    m.addConstr(sum(x[i,j] for j in suppliers) == demand[i])
for j in suppliers:
    m.addConstr(sum(x[i,j] for i in customers) <= capacity[j]*y[j])
    
# m.Params.OutputFlag = 0
m.optimize()
m.Params.OutputFlag = 1
m.printAttr('X')

Set parameter Username
Academic license - for non-commercial use only - expires 2023-11-04
Gurobi Optimizer version 9.5.2 build v9.5.2rc0 (win64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 12 rows, 36 columns and 68 nonzeros
Model fingerprint: 0x5385f77b
Variable types: 32 continuous, 4 integer (4 binary)
Coefficient statistics:
  Matrix range     [1e+00, 7e+02]
  Objective range  [1e+01, 5e+03]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+02, 2e+02]
Presolve time: 0.01s
Presolved: 12 rows, 36 columns, 68 nonzeros
Variable types: 32 continuous, 4 integer (4 binary)
Found heuristic solution: objective 34925.000000

Root relaxation: objective 2.618097e+04, 7 iterations, 0.00 seconds (0.00 work units)

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0 26180.9740    0    4 34925.0000 26180.9740  25.0%     -    0s
H    0     0                    29280.000000 26180.9740  10.6%     -    0s
     0     0 28490.0000    0    3 29280.0000 28490.0000  2.70%     -    0s
     0     0 29280.0000    0    2 29280.0000 29280.0000  0.00%     -    0s

Cutting planes:
  Gomory: 2
  Implied bound: 4

Explored 1 nodes (19 simplex iterations) in 0.08 seconds (0.00 work units)
Thread count was 16 (of 16 available processors)

Solution count 2: 29280 34925 

Optimal solution found (tolerance 1.00e-04)
Best objective 2.928000000000e+04, best bound 2.928000000000e+04, gap 0.0000%

    Variable            X 
-------------------------
       y[S1]            1 
       y[S2]            1 
    x[D1,S1]           10 
    x[D1,S2]           90 
    x[D2,S1]          150 
    x[D3,S1]          175 
    x[D4,S1]          125 
    x[D5,S2]          180 
    x[D6,S1]          140 
    x[D7,S2]          120 
    x[D8,S2]          160 

The lazy constraint callback will add lazy cut for any MIP solution found during the branch-and-bound tree. Hence, the list of objective value is not necessary monotonic.

Code

master = Model()
y = master.addVars(suppliers, vtype=GRB.BINARY, name='y')
z = master.addVar(vtype=GRB.CONTINUOUS, name='z')
master.setObjective(sum(fixed_cost[j]*y[j] for j in suppliers) + z)
# master.addConstr(sum(capacity[j]*y[j] for j in suppliers) >= sum(demand[i] for i in customers))
# _vars is used in callback function
master._vars = {'y': y, 'z': z}

# primal formulation of the subproblem
subP = Model()
x = subP.addVars(customers, suppliers, vtype=GRB.CONTINUOUS, name='x')
subP.setObjective(sum(trans_cost[i,j]*x[i,j] for i in customers for j in suppliers))
for i in customers:
    subP.addConstr(sum(x[i,j] for j in suppliers) == demand[i], f"d[{i}]")
for j in suppliers:
    subP.addConstr(sum(x[i,j] for i in customers) <= capacity[j]*1, f"c[{j}]")
subP.Params.InfUnbdInfo = 1
# have to update (otherwise empty) so that we can get constr and its rhs later in
#   subP.getConstrByName(f"c[{j}]").RHS
subP.update()

# dual formulation of the subproblem
subD = Model()
u = subD.addVars(customers, lb=-1e20, ub=1e20, vtype=GRB.CONTINUOUS, name='u')
v = subD.addVars(suppliers, vtype=GRB.CONTINUOUS, name='v')
for i in customers:
    for j in suppliers:
        subD.addConstr(u[i] - v[j] <= trans_cost[i,j])
subD.Params.InfUnbdInfo = 1

master.Params.OutputFlag = 0
subP.Params.OutputFlag = 0
subD.Params.OutputFlag = 0

def benders():
    objs = []
    niter = 0
    _USE_SUB_PRIMAL = True
    while True:
        master.optimize()
        objs.append(master.objVal)

        if _USE_SUB_PRIMAL:
            for j in suppliers:
                subP.getConstrByName(f"c[{j}]").RHS = capacity[j]*y[j].x
            subP.optimize()

            if subP.status == GRB.INFEASIBLE:
                u_ray = dict([(i, -subP.getConstrByName(f"d[{i}]").FarkasDual) for i in customers])
                v_ray = dict([(j, subP.getConstrByName(f"c[{j}]").FarkasDual) for j in suppliers])

                exp = LinExpr(sum(demand[i]*u_ray[i] for i in customers) -
                              sum(capacity[j]*y[j]*v_ray[j] for j in suppliers))
                master.addConstr(exp <= 0) 

            if subP.status == GRB.OPTIMAL:
                u_ext = dict([(i, subP.getConstrByName(f"d[{i}]").Pi) for i in customers])
                v_ext = dict([(j, -subP.getConstrByName(f"c[{j}]").Pi) for j in suppliers])

                exp = LinExpr(sum(demand[i]*u_ext[i] for i in customers) -
                              sum(capacity[j]*y[j]*v_ext[j] for j in suppliers))
                if exp.getValue() > z.x:
                    master.addConstr(exp <= z) 
                else:
                    break    
        else:
            subD.setObjective(sum(demand[i]*u[i] for i in customers) - 
                              sum(capacity[j]*y[j].x*v[j] for j in suppliers), GRB.MAXIMIZE)
            subD.optimize()

            if subD.status == GRB.UNBOUNDED:
                master.addConstr(sum(demand[i]*u[i].unbdray for i in customers) - 
                                 sum(capacity[j]*y[j]*v[j].unbdray for j in suppliers) <= 0)

            if subD.status == GRB.OPTIMAL:
                exp = LinExpr(sum(demand[i]*u[i].x for i in customers) - 
                              sum(capacity[j]*y[j]*v[j].x for j in suppliers))
                if exp.getValue() > z.x:
                    master.addConstr(exp <= z) 
                else:
                    break        
    #     if subP.status == GRB.INFEASIBLE:            
    #         print([u[i].unbdray for i in customers])
    #         print([v[j].unbdray for j in suppliers])
    #         print([subP.getConstrByName(f"d[{i}]").FarkasDual for i in customers])
    #         print([subP.getConstrByName(f"c[{j}]").FarkasDual for j in suppliers])
    #     if subP.status == GRB.OPTIMAL:
    #         print([u[i].x for i in customers])
    #         print([v[j].x for j in suppliers])
    #         print([subP.getConstrByName(f"d[{i}]").Pi for i in customers])
    #         print([subP.getConstrByName(f"c[{j}]").Pi for j in suppliers])
        niter = niter + 1
    return objs

niter = 0
objs = []
def cbbenders(model, where):
    global niter
    if where == GRB.Callback.MIPSOL:
        # print('-'*8 + " MIP sol callback: {} ".format(niter) + '-'*8)

        obj = model.cbGet(GRB.Callback.MIPSOL_OBJ)
        objs.append(round(obj))

        y_val = model.cbGetSolution(model._vars['y'])
        z_val = model.cbGetSolution(model._vars['z'])
        subD.setObjective(sum(demand[i]*u[i] for i in customers) -
                          sum(capacity[j]*y_val[j]*v[j] for j in suppliers), GRB.MAXIMIZE)
        # turn on option to query unbounded ray
        subD.Params.InfUnbdInfo = 1
        subD.update()
        subD.optimize()
        
        if subD.status == GRB.UNBOUNDED:
            master.cbLazy(sum(demand[i]*u[i].unbdray for i in customers) - 
                             sum(capacity[j]*model._vars['y'][j]*v[j].unbdray for j in suppliers) <= 0)

        if subD.status == GRB.OPTIMAL:
            rhs = (sum(demand[i]*u[i].x for i in customers) - 
                   sum(capacity[j]*y_val[j]*v[j].x for j in suppliers))
            if rhs > z_val:
                master.cbLazy(sum(demand[i]*u[i].x for i in customers) - 
                              sum(capacity[j]*model._vars['y'][j]*v[j].x for j in suppliers) <= model._vars['z']) 
        niter = niter+1
        # print('-'*8 + " MIP sol callback: End " + '-'*8)

master.Params.lazyConstraints = 1 
# master.Params.PreCrush = 1
# master.Params.Threads = 1

# Direct Implementation
# benders()

# Lazy Constraint
master.Params.OutputFlag = 1
master.optimize(cbbenders)
print(objs)

Set parameter InfUnbdInfo to value 1
Set parameter InfUnbdInfo to value 1
Set parameter OutputFlag to value 1
Gurobi Optimizer version 9.5.2 build v9.5.2rc0 (win64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 0 rows, 5 columns and 0 nonzeros
Model fingerprint: 0xd1b1e689
Variable types: 1 continuous, 4 integer (4 binary)
Coefficient statistics:
  Matrix range     [0e+00, 0e+00]
  Objective range  [1e+00, 5e+03]
  Bounds range     [1e+00, 1e+00]
  RHS range        [0e+00, 0e+00]
Presolve time: 0.00s
Presolved: 0 rows, 5 columns, 0 nonzeros
Variable types: 1 continuous, 4 integer (4 binary)

Root relaxation: objective 7.861538e+03, 1 iterations, 0.00 seconds (0.00 work units)

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0 7861.53846    0    1          - 7861.53846      -     -    0s
H    0     0                    30560.000000 7861.53846  74.3%     -    0s
H    0     0                    29280.000000 7861.53846  73.2%     -    0s
     0     0 22620.2509    0    2 29280.0000 22620.2509  22.7%     -    0s
     0     0 25361.1345    0    4 29280.0000 25361.1345  13.4%     -    0s
     0     0 25687.3911    0    4 29280.0000 25687.3911  12.3%     -    0s
     0     0 25717.6280    0    4 29280.0000 25717.6280  12.2%     -    0s
     0     0 25890.4329    0    4 29280.0000 25890.4329  11.6%     -    0s
     0     0 25908.8205    0    4 29280.0000 25908.8205  11.5%     -    0s
     0     0 25993.6985    0    4 29280.0000 25993.6985  11.2%     -    0s
     0     0 25995.8961    0    4 29280.0000 25995.8961  11.2%     -    0s
     0     0 26086.8822    0    4 29280.0000 26086.8822  10.9%     -    0s
     0     0 26219.4737    0    4 29280.0000 26219.4737  10.5%     -    0s
     0     0 26222.6884    0    4 29280.0000 26222.6884  10.4%     -    0s
     0     0 26334.1892    0    4 29280.0000 26334.1892  10.1%     -    0s
     0     0 26335.7739    0    4 29280.0000 26335.7739  10.1%     -    0s
     0     0 26375.6269    0    4 29280.0000 26375.6269  9.92%     -    0s
     0     0 26377.0966    0    4 29280.0000 26377.0966  9.91%     -    0s
     0     0 26377.3105    0    4 29280.0000 26377.3105  9.91%     -    0s
     0     0 26396.4657    0    4 29280.0000 26396.4657  9.85%     -    0s
     0     0 26400.7512    0    4 29280.0000 26400.7512  9.83%     -    0s
     0     2 26400.7512    0    4 29280.0000 26400.7512  9.83%     -    0s

Cutting planes:
  MIR: 4
  Lazy constraints: 7

Explored 9 nodes (49 simplex iterations) in 0.03 seconds (0.00 work units)
Thread count was 16 (of 16 available processors)

Solution count 2: 29280 30560 

Optimal solution found (tolerance 1.00e-04)
Best objective 2.928000000000e+04, best bound 2.928000000000e+04, gap 0.0000%

User-callback calls 186, time in user-callback 0.01 sec
[0, 0, 0, 0, 0, 8650, 30560, 24210, 29280, 26251, 25630, 28460, 27860]

Dantzig-Wolfe Decomposition

Multicommodity Flow Problem

\begin{align} &\text{minimize} &\sum_{k \in K} \sum_{(i,j) \in A} c_{ij}^k x_{ij}^k \nonumber \\ &\text{subject to} &\sum_{k \in K} x_{ij}^k &\le u_{ij} && (i,j) \in A \nonumber \\ &&\sum_{(i,j) \in A} x_{ij}^k - \sum_{(j,i) \in A} x_{ji}^k &= b_i^k && i \in N,\ k \in K \\ && 0 \le x_{ij}^k &\le u_{ij}^k && (i,j) \in A,\ k \in K \end{align}

Matrix Form

The matrix form of the constraints is:

\begin{equation*} \scriptsize \left( \begin{array} {rrrrrrr|rrrrrrr} x^1_{12} & & & & & & & x^2_{12} & & & & & \\ & x^1_{13} & & & & & & & x^2_{13} & & & & \\ & & x^1_{53} & & & & & & & x^2_{53} & & & \\ & & & x^1_{56} & & & & & & & x^2_{56} & & \\ & & & & x^1_{34} & & & & & & & x^2_{34} & \\ & & & & & x^1_{24} & & & & & & & x^2_{24}& \\ & & & & & & x^1_{46} & & & & & & & x^2_{46}\\ \hline x^1_{12} & x^1_{13} & & & & & \\ -x^1_{12} & & & & & -x^1_{24} & \\ & -x^1_{13} & -x^1_{53} & & x^1_{34} & & \\ & & & & -x^1_{34} & x^1_{24} & x^1_{46} \\ & & x^1_{53} & x^1_{56} & & & \\ & & & -x^1_{56} & & & -x^1_{46} \\ \hline && & & & & & x^2_{12} & x^2_{13} & & & & & & \\ && & & & & & -x^2_{12} & & & & & -x^2_{24} & & \\ && & & & & & & -x^2_{13} & -x^2_{53} & & x^2_{34} & & \\ && & & & & & & & & & -x^2_{34} & x^2_{24} & x^2_{46} \\ && & & & & & & & x^2_{53} & x^2_{56} & & & & \\ && & & & & & & & & -x^2_{56} & & & -x^2_{46} & \end{array} \right) \end{equation*}

The coefficient matrix is \begin{equation*} \left( \begin{array} {rrrrrrr|rrrrrrr} 1 & & & & & & & 1 \\ & 1 & & & & & & & 1 \\ & & 1 & & & & & & & 1 \\ & & & 1 & & & & & & & 1 \\ & & & & 1 & & & & & & & 1 \\ & & & & & 1 & & & & & & & 1 \\ & & & & & & 1 & & & & & & & 1 \\ \hline 1 & 1 & & & & & \\ -1 & & & & & -1 & \\ & -1 & -1 & & 1 & & \\ & & & & -1 & 1 & 1 \\ & & 1 & 1 & & & \\ & & & -1 & & & -1 \\ \hline & & & & & & & 1 & 1 \\ & & & & & & & -1 & & & & & -1 \\ & & & & & & & & -1 & -1 & & 1 \\ & & & & & & & & & & & -1 & 1 & 1 \\ & & & & & & & & & 1 & 1 \\ & & & & & & & & & & -1 & & & -1 \end{array} \right) \end{equation*}

\begin{equation*} \begin{pmatrix} A_1 & A_2 & \cdots & A_{|K|} \\ \hline B_1 & & & \\ & B_2 & & \\ & & \ddots & \\ & & & B_{|K|} \end{pmatrix} \end{equation*}

The Multicommodity Flow Problem can be formulated as follows: \begin{align*} &\text{minimize} &\sum_{k \in K} \mathbf{c}^k \mathbf{x}^k \nonumber \\ &\text{subject to} &\sum_{k \in K} (\mathbf{x}^k)_{ij} &\le u_{ij} && (i,j) \in A \\ &&\mathbf{x}^k &\in X^k && k \in K \end{align*}

where X^k = \{\mathbf{x}^k \in \mathbb{R}^{|A|}: \text{$\mathbf{x}^k$ satisfies (1)-(2)}\}. Let \left\{\mathbf{x}_p^k\right\}_{p \in P^k} be all extreme points of X^k (where p stands for a proposal). We have \begin{align*} \mathbf{x}^k \in X^k \Leftrightarrow \mathbf{x}^k = \sum\limits_{p \in P^k} \lambda_p^k \mathbf{x}_p^k, \sum\limits_{p \in P^k} \lambda_p^k = 1, \lambda_p^k \ge 0 \end{align*}

Similarly, - c_p^k = \mathbf{c}^k \mathbf{x}_p^k is cost of proposal p \in P^k - f_{pij}^k = (\mathbf{x}_p^k)_{ij} is flow of commodity k along (i,j) for p \in P^k

The Multicommodity Flow Problem can be formulated as follows: \begin{align*} \text{minimize} &\sum_{k \in K} \sum_{p \in P^k} c_p^k \lambda_p^k \nonumber \\ \text{subject to} &\sum_{k \in K} \sum_{p \in P^k} f^k_{pij} \lambda_p^k \le u_{ij} && (i,j) \in A && (v_{ij} \le 0)\\ &\sum_{p \in P^k} \lambda_p^k = 1 && k \in K && (\text{$w_k$ free})\\ &\lambda_p^k \ge 0 && k \in K,\ p \in P^k \end{align*}

Code

# Base data
commodities = [1, 2]
nodes = [1, 2, 3, 4, 5, 6]

arcs, capacity = multidict({
    (1, 2): 5,
    (1, 3): 30,
    (3, 4): 10,
    (4, 2): 30,
    (4, 6): 30,
    (5, 3): 30,
    (5, 6): 30})

# Cost for triplets commodity-source-destination
cost = {
    (1, 1, 2): 1,
    (1, 1, 3): 5,
    (1, 5, 3): 1,
    (1, 5, 6): 5,
    (1, 3, 4): 1,
    (1, 4, 2): 5,
    (1, 4, 6): 1,
    (2, 1, 2): 1,
    (2, 1, 3): 5,
    (2, 5, 3): 1,
    (2, 5, 6): 5,
    (2, 3, 4): 1,
    (2, 4, 2): 5,
    (2, 4, 6): 1}

# Supply for pairs of commodity-node
outflow = {
    (1, 1):  10,
    (1, 2): -10,
    (1, 3):  0,
    (1, 4):  0,
    (1, 5):  0,
    (1, 6):  0,
    (2, 1):  0,
    (2, 2):  0,
    (2, 3):  0,
    (2, 4):  0,
    (2, 5):  20,
    (2, 6): -20}

m = Model()
x = m.addVars(commodities, arcs, obj=cost, name="flow")

# Arc-capacity constraints
m.addConstrs(
    (x.sum('*', i, j) <= capacity[i, j] for i, j in arcs), "cap")

# Flow-conservation constraints
m.addConstrs(
    (x.sum(k, i, '*') - x.sum(k, '*', i) == outflow[k, i]
        for k in commodities for i in nodes), "node")
            
m.optimize()
# m.write('test.lp')

Gurobi Optimizer version 9.5.2 build v9.5.2rc0 (win64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 19 rows, 14 columns and 42 nonzeros
Model fingerprint: 0xa507224a
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  Objective range  [1e+00, 5e+00]
  Bounds range     [0e+00, 0e+00]
  RHS range        [5e+00, 3e+01]
Presolve removed 19 rows and 14 columns
Presolve time: 0.00s
Presolve: All rows and columns removed
Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    1.5000000e+02   0.000000e+00   0.000000e+00      0s

Solved in 0 iterations and 0.00 seconds (0.00 work units)
Optimal objective  1.500000000e+02

Code

master = Model()
excess = master.addVar(0.0, GRB.INFINITY, 0.0, GRB.CONTINUOUS)
master.addConstrs((-excess <= capacity[i,j] for (i,j) in arcs), "multi")
master.addConstr(0*excess == 1, "convex")

sub = Model()
x = sub.addVars(commodities, arcs, vtype=GRB.CONTINUOUS, name="flow")
sub.addConstrs((x.sum(k, i, '*') - x.sum(k, '*', i) == outflow[k, i] for k in commodities for i in nodes), "node")

master.setObjective(excess, GRB.MINIMIZE)
sub.setObjective(quicksum(0*x[k,i,j]-1 for k in commodities for (i,j) in arcs), GRB.MINIMIZE)

nprop = 0

prop_cost = dict({})
prop_ship = dict({})

master.Params.OutputFlag = 0
sub.Params.OutputFlag = 0

while True:
    sub.optimize()

    if sub.objVal >= -1e-5:
        print('NO FEABIEL SOLUTION')
        break
    
    prop_cost[nprop] = sum(cost[k,i,j]*x[k,i,j].x for k in commodities for (i,j) in arcs)
    prop_ship[nprop] = dict([((i,j),sum(x[k,i,j].x for k in commodities)) for (i,j) in arcs])
    print(prop_ship)
    master.update() # must update otherwise master.getConstrs() gets empty
    master.addVar(vtype=GRB.CONTINUOUS, name='l['+str(nprop)+']',
                  column=Column([*prop_ship[nprop].values(),1], master.getConstrs()))
    master.optimize()
    nprop = nprop+1
    if master.objVal <= 1e-5:
        break
    else:
        v = dict([((i,j),master.getConstrByName(f"multi[{i},{j}]").Pi) for (i,j) in arcs])
        w = master.getConstrByName(f"convex").Pi
        sub.setObjective(quicksum((-v[i,j])*x[k,i,j] for k in commodities for (i,j) in arcs) - 
                         sum(1 for k in commodities), GRB.MINIMIZE)

print(v,w)        
master.Params.OutputFlag = 1
sub.Params.OutputFlag = 1

excess.lb = 0
excess.ub = 0
master.setObjective(quicksum(prop_cost[t]*master.getVarByName(f"l[{t}]") for t in range(nprop)), GRB.MINIMIZE)
master.update()

while True:
    sub.setObjective(quicksum((cost[k,i,j]-v[i,j])*x[k,i,j] for k in commodities for (i,j) in arcs) -
                     sum(w for k in commodities), GRB.MINIMIZE)
    sub.optimize()

    if sub.objVal >= -1e-5:
        print('OPTIMAL SOLUTION')
    else:
        break

    prop_cost[nprop] = sum(cost[k,i,j]*x[k,i,j].x for k in commodities for (i,j) in arcs)
    prop_ship[nprop] = dict([((i,j),sum(x[k,i,j].x for k in commodities)) for (i,j) in arcs])
    master.update() # must update otherwise master.getConstrs() gets empty
    master.addVar(vtype=GRB.CONTINUOUS,  name='l['+str(nprop)+']', obj=prop_cost[nprop],
                  column=Column([*prop_ship[nprop].values(),1], master.getConstrs()))
    master.optimize()
    nprop = nprop+1
    v = dict([((i,j),master.getConstrByName(f"multi[{i},{j}]").Pi) for (i,j) in arcs])
    w = master.getConstrByName(f"convex").Pi
    
opt_ship = dict([((i,j), sum(prop_ship[t][i,j]*master.getVarByName(f"l[{t}]").x for t in range(nprop))) for (i,j) in arcs])
# opt_ship

{0: {(1, 2): 10.0, (1, 3): 0.0, (5, 3): 0.0, (5, 6): 20.0, (3, 4): 0.0, (4, 2): 0.0, (4, 6): 0.0}}
{0: {(1, 2): 10.0, (1, 3): 0.0, (5, 3): 0.0, (5, 6): 20.0, (3, 4): 0.0, (4, 2): 0.0, (4, 6): 0.0}, 1: {(1, 2): 0.0, (1, 3): 10.0, (5, 3): 0.0, (5, 6): 20.0, (3, 4): 10.0, (4, 2): 10.0, (4, 6): 0.0}}
{(1, 2): -1.0, (1, 3): 0.0, (5, 3): 0.0, (5, 6): 0.0, (3, 4): 0.0, (4, 2): 0.0, (4, 6): 0.0} 10.0
Set parameter OutputFlag to value 1
Set parameter OutputFlag to value 1
Gurobi Optimizer version 9.5.2 build v9.5.2rc0 (win64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 12 rows, 14 columns and 28 nonzeros
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  Objective range  [1e+00, 5e+00]
  Bounds range     [0e+00, 0e+00]
  RHS range        [1e+01, 2e+01]
Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0   -1.1000000e+31   4.000000e+30   1.100000e+01      0s
       2    6.0000000e+01   0.000000e+00   0.000000e+00      0s

Solved in 2 iterations and 0.00 seconds (0.00 work units)
Optimal objective  6.000000000e+01
OPTIMAL SOLUTION
Gurobi Optimizer version 9.5.2 build v9.5.2rc0 (win64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 8 rows, 4 columns and 20 nonzeros
Coefficient statistics:
  Matrix range     [1e+00, 2e+01]
  Objective range  [7e+01, 2e+02]
  Bounds range     [0e+00, 0e+00]
  RHS range        [1e+00, 3e+01]
Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0   -6.0000000e+31   4.500000e+30   6.000000e+01      0s
       3    1.5000000e+02   0.000000e+00   0.000000e+00      0s

Solved in 3 iterations and 0.00 seconds (0.00 work units)
Optimal objective  1.500000000e+02
Gurobi Optimizer version 9.5.2 build v9.5.2rc0 (win64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 12 rows, 14 columns and 28 nonzeros
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  Objective range  [1e+00, 1e+01]
  Bounds range     [0e+00, 0e+00]
  RHS range        [1e+01, 2e+01]
Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0   -2.3000000e+02   0.000000e+00   0.000000e+00      0s

Solved in 0 iterations and 0.00 seconds (0.00 work units)
Optimal objective -2.300000000e+02

Code

m = Model()
x = m.addVars(commodities, arcs, vtype=GRB.CONTINUOUS, name="flow")
m.setObjective(quicksum(cost[k,i,j]*x[k,i,j] for k in commodities for (i,j) in arcs), GRB.MINIMIZE)
m.addConstrs((x.sum(k, i, '*') - x.sum(k, '*', i) == outflow[k, i] for k in commodities for i in nodes), "node")
m.addConstrs(quicksum(x[k,i,j] for k in commodities) == opt_ship[i,j] for (i,j) in arcs)
m.optimize()
x

Gurobi Optimizer version 9.5.2 build v9.5.2rc0 (win64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 19 rows, 14 columns and 42 nonzeros
Model fingerprint: 0xd534693b
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  Objective range  [1e+00, 5e+00]
  Bounds range     [0e+00, 0e+00]
  RHS range        [5e+00, 2e+01]
Presolve removed 19 rows and 14 columns
Presolve time: 0.00s
Presolve: All rows and columns removed
Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    1.5000000e+02   0.000000e+00   0.000000e+00      0s

Solved in 0 iterations and 0.00 seconds (0.00 work units)
Optimal objective  1.500000000e+02

{(1, 1, 2): <gurobi.Var flow[1,1,2] (value 5.0)>,
 (1, 1, 3): <gurobi.Var flow[1,1,3] (value 5.0)>,
 (1, 5, 3): <gurobi.Var flow[1,5,3] (value 0.0)>,
 (1, 5, 6): <gurobi.Var flow[1,5,6] (value 0.0)>,
 (1, 3, 4): <gurobi.Var flow[1,3,4] (value 5.0)>,
 (1, 4, 2): <gurobi.Var flow[1,4,2] (value 5.0)>,
 (1, 4, 6): <gurobi.Var flow[1,4,6] (value 0.0)>,
 (2, 1, 2): <gurobi.Var flow[2,1,2] (value 0.0)>,
 (2, 1, 3): <gurobi.Var flow[2,1,3] (value 0.0)>,
 (2, 5, 3): <gurobi.Var flow[2,5,3] (value 5.0)>,
 (2, 5, 6): <gurobi.Var flow[2,5,6] (value 15.0)>,
 (2, 3, 4): <gurobi.Var flow[2,3,4] (value 5.0)>,
 (2, 4, 2): <gurobi.Var flow[2,4,2] (value 0.0)>,
 (2, 4, 6): <gurobi.Var flow[2,4,6] (value 5.0)>}

Code

# sub.display()
# sub.getVars()
master.display()
master.getVars()

Minimize
  <gurobi.LinExpr: 110.0 l[0] + 210.0 l[1] + 70.0 l[2]>
Subject To
  multi[1,2]: <gurobi.LinExpr: -1.0 C0 + 10.0 l[0] + 10.0 l[2]> <= 5
  multi[1,3]: <gurobi.LinExpr: -1.0 C0 + 10.0 l[1]> <= 30
  multi[5,3]: <gurobi.LinExpr: -1.0 C0 + 20.0 l[2]> <= 30
  multi[5,6]: <gurobi.LinExpr: -1.0 C0 + 20.0 l[0] + 20.0 l[1]> <= 30
  multi[3,4]: <gurobi.LinExpr: -1.0 C0 + 10.0 l[1] + 20.0 l[2]> <= 10
  multi[4,2]: <gurobi.LinExpr: -1.0 C0 + 10.0 l[1]> <= 30
  multi[4,6]: <gurobi.LinExpr: -1.0 C0 + 20.0 l[2]> <= 30
  convex: <gurobi.LinExpr: l[0] + l[1] + l[2]> = 1
Bounds
  C0 = 0

[<gurobi.Var C0 (value 0.0)>,
 <gurobi.Var l[0] (value 0.25)>,
 <gurobi.Var l[1] (value 0.5)>,
 <gurobi.Var l[2] (value 0.25)>]

GAP

\begin{align*} &\text{maximize} &\sum_{1 \le i \le m} \sum_{1 \le j \le n} f_{ij} x_{ij} \\ &\text{subject to} &\sum_{1 \le i \le m} x_{ij} & = 1 && j = 1 \ldots n \\ &&\sum_{1 \le j \le n} w_{ij} x_{ij} &\le c_i && i = 1 \ldots m \\ && x_{ij} & \in \{0,1\} && i = 1 \ldots m; j = 1 \ldots n \end{align*}

• x^i_p = (x^i_{1p},\ldots, x^i_{np}): a solution to agent i such that

  • x^i_{jp}=1: agent i is assigned to job j in the pth (this) solution

  • x^i_{jp}=0 otherwise;

• P_i = \{x^i_1, \ldots, x^i_{p_i}\}: set of all possible feasible assignments for agent i;

• \lambda^i_p \in \color{red}{\{0,1\}}: whether a feasible assignment x^i_p is selected;

• x_{ij} = \sum_{1 \le p \le p_i} x^i_{jp} \lambda^i_p.

\begin{align*} &\text{maximize} &\sum_{1 \le i \le m} \sum_{1 \le p \le p_i} \left( \sum_{1 \le j \le n} f_{ij} x^i_{jp} \right) \lambda^i_p\\ &\text{subject to} &\sum_{1 \le i \le m} \sum_{1 \le p \le p_i} x^i_{jp} \lambda^i_p & = 1 && j = 1 \ldots n \\ &&\sum_{1 \le p \le p_i} \lambda^i_{p} &\le 1 && i = 1 \ldots m \\ && \lambda^i_{p} & \in \{0,1\} && i = 1 \ldots m; p \in P_i \end{align*}